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Introduction 


Work under Cooperative Agreement NCC2-312 was performed 
in two stages. Stage one commenced on 1 June 1984 and 
extended through 2/28/87. It was performed under the 
direction of Dr. Bala A. Balakrishnan as Eloret Principal 
Investigator; Mr. Mike Green of NASA-Ames Research Center 
being the Technical Monitor. The final reporting for that 
stage of the research program was submitted on 31 July, 1987. 
It included the following five AIAA papers: 85-1006; 85-1064; 

85-1063; 86-1277; and 86-1312. 

Stage two was performed under the Direction of Dr. Dinesh 
K. Prabhu during the period 7/1/87 through 5/31/90. Mr. 
William C. Davy and Dr. George S. Deiwert were the NASA 
Technical Monitors. This final report concentrates on Dr. 
Prabhu ' s work. 
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Abstract 


The equations governing the multidimensional flow of a reacting 
mixture of thermally perfect gases have been derived. The modelling 
procedures for the various terms of the conservation laws have been 
discussed. A numerical algorithm (based on the finite-volume approach) 
to solve these conservation equations has been developed. The ad- 
vantages and disadvantages of the present numerical scheme have been 
discussed from the point of view of accuracy, computer time and memory 
requirements. A simple one-dimensional model problem has been solved 
to prove the feasibility and accuracy fo the algorithm. A computer 
code implementing the above algorithm has been developed and is cur- 
rently being applied to simple geometries and comditions. Once the 
code is completely debugged and validated it will be used to compute 
the complete, unsteady flow field aournd the Aeroassist Flight Exper- 
iment (AFE) body. 

Nomenclature 


0 : null vector 

aj : frozen speed of sound 

c a : mass fraction of species s 

C V j : frozen specific heat of the mixture 

C v a : specific heat of species s 

D : binary diffusion coefficient 

E • total energy per unit volume of the mixture 

e : specific internal energy of the mixture 

e s : specific internal energy of species s 
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: algebraic flux vector 

: numerical flux vector 

: normal component of the flux 

: specific total enthalpy of the mixture 

; specific enthalpy of formation of species s 
identity tensor 

: backward reaction rate constant for the 771th reaction 

: forward reaction rate constant for the 771th reaction 

: reference length (777) 

: binary Lewis number 

: number of reactions 

: frozen Mach number 

: molecular weight 

number of dimensions 
: number of conservation equations 

: number of species 

unit normal to cell surface 
: static pressure 

: pressure tensor 

: heat flux vector 

: algebraic vector of conservative variables 

: universal gas constant, 8314.34 J/ikmol-K') 

: position vector 

: Reynolds number based on L 

: magnitude of surface area vector 

: unit tangent to cell surface 
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S : 

surface area vector 

t : 

time 

t 

unit tangent to cell surface 

T : 

temperature 

T : 

viscous stress tensor 

U 

algebraic vector of primitive variables 

u : 

mass-averaged velocity 

«n : 

velocity component normal to cell surface 

V : 

magnitude of the mass— averaged velocity 

V : 

cell volume 

V, : 

diffusion velocity of species s 

v : 

control surface velocity 

w a : 

mass production rate of species s 

W : 

vector of production terms 

W : 

numerical representation of the source vector 

: 

mole fraction of species s 


molar concentration per unit volume of species s 

K • 

thermal conductivity of the mixture 

K a : 

thermal conductivity of species s 

/* : 

viscosity of the mixture 

Ha : 

viscosity of species s 

P : 

mass density of the mixture 

P* : 

mass density of species s 

Subscripts 


j,k,l ■■ 

cell centroid indices 

y,i,o±l : 

cell face indices 
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n 

r,s,l 
ref 
w 
oo 

Superscripts 


* 

v 


a 


time index 

indices denoting species 
reference condition 
wall 

f reestream 

dimensional quantity 
inviscid/convect ive quantity 
viscous /nonconvect ive quantity 
direction index 


Governing Equations 

Let V be the control volume of interest and dV its bounding control 
surface which moves with a velocity v . The equations governing the 
flow are 1 ’ 2 : 

(a) the law of conservation of mass of species s (s = 1, 2, . . . ,n 3 ) : 

— [ p.dV + { pJu. - v) • ndA + l p s V a -ndA= w s d.V (1) 

dtJv(t) JdV(t) JdV(t) Mt) 

The first term on the left hand side of Eq. 1 represents the time-rate 

of change of species mass in the control volume, the second term rep 

resents the net flux across the control surface and the third term 

represents the diffusion mass flux. The term on the right hand side 

of Eq. 1 represents the rate of mass production/depletion within the 

control volume. The equation for the overall conservation of mass is 

obtained by summing Eq. 1 over all the species and noting that 

£>.=* E^ v *=°- I>* = ° 

. 8 a 
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Thus, the global continuity equation is 


— [ pdV+ l p { u - v) • ndA = 0 

diJv(t) JdV(t) 


I V(t) JdV(t) 

(b) the law of conservation of linear momentum: 


— f pndV + i pu(u - v) • ndA - <f P ndA = 0 
dtjv(t) JdVit) JdV(t) 


( 2 ) 


/V(t) JdV{t) JdV(t) 

where body forces have been neglected and P represents the surface 
forces (stress tensor) . The above vector equation leads to n d scalar 
equations, where n d is the number of dimensions. 

(c) the law of conservation of energy: 


— [ EdV + i E(u - v) • ndA = ~l q * ndA + i P • u • ndA 
dt Jv(t) Lv(t) dav(t) Jam 

where E is volumetric total energy of the mixture, and q is the heat 
flux vector. The second term on the right hand side of Eq. 3 repre 
sents the work done by the surface forces. The volumetric total en- 
ergy of the mixture is defined as 


( 3 ) 


E = p(e + -u-u) 


( 4 ) 


We have, therefore, a total of n e = n s + n d + 1 conservation equations 
to solve. The scalar and vector fields of interest to us are p a , p, u, 
p, and T. In order to close the set of governing equations we need to 
model the various terms appearing in the set. This is discussed in the 
following section. 
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Physical Model 


(a) thermal equation of state : 

The constituent gases of the mixture are assumed, to be thermally 
perfect and obey the following equation of state for the partial pres- 
sure p a 

p s 1lT 

Pa = ^l 7 

Using Dalton' s law of partial pressures, the following equation of 
state is obtained for the mixture 


P = 


pTZT 

M 


(5) 


where M is the molecular weight of the mixture and is defined as 


*■(?*) 


-1 


(7) 


(b) thermodynamic properties : 

The thermodynamic properties of the individual species of the gas 
are available as cubic-spline fits. These curve fits are of the form 


e a = e s (T ) + h° 3 = ao + a\T + a 2 T 2 + a 3 T 3 
C„. = C„.(T) = ^ = a, + 2o 2 T + 3a 3 T 2 

where a 0 ,...,a 3 are the spline coefficients. These coefficients have 
been generated by Liu and Vinokur 3 and based on quantum calculations 
done by Jaf fe 4 . These thermodynamic properties are considered more 


7 



realistic than previous calculations because they account for the in 
ternal structure of the atoms and molecules. The thermodynamic prop- 
erties of mixture as a whole are computed as simple weighted sums of 
the individual species properties, i.e., 

C Vj =Y, CsCv ‘ e = £c.e. (8 

S 3 

(C) transport properties : 

The viscosities of the individual species are computed using curve 
fits developed by Blottner et al . 5 and these are of the form 


= exp[(a log e T + b) log e T + c] 


where a,b,c are constants. The thermal conductivities of the individ 
ual species are computed using Eucken's formula 

3 m s k n 


The transport properties of the mixture are calculated using Wilke s 
mixing rule . 6 This mixing rule, considered adequate for weakly ioniz- 
ing flows, is mathematically expressed as 


f-E 



-E 



( 9 ) 


where 


X 3 = 


c 3 M 

M a 


*. = 


r— 1 


1 + 
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The binary diffusion coefficient is obtained from the definition 
of the Lewis number, £e, which is assumed to be the same for all species . 

Therefore, 


V = 


n£e 

pC Pj 


( 10 ) 


(d) diffusion mass flux : 

Assuming the thermal and pressure diffusion effects to be negligi- 
ble and that mass diffusion is binary, the diffusion mass flux can be 
expressed as 

P.V, = -pVVc, (11) 


The binary diffusion approximation is valid for mildly ionizing air 
consisting of "heavy" particles (molecules and ions) and "light" par- 
ticles (atoms) . This assumption considerably simplifies the complex- 
ities of a true multicomponent diffusion model. 

(e) stress tensor : 

Bulk viscosity effects are considered to be negligibly small . The 
stress tensor represented as the sum of the hydrostatic and devia- 
toric stress is written as 

P = —pi + T = —pi + |i[Vu + (Vu) T - |v ■ ul] (12) 


where p is the coefficient of viscosity, I is the identity tensor, and 
V is the usual gradient operator . 

(f) heat flux vector : 

The heat flux vector consists of two parts, (i) a conduction part 
and (ii) a diffusion part. This is expressed as 


q = — kVT + 'y haPsV a 

3 


( 13 ) 
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where k is the coefficient of thermal conductivity. The heat transfer 
due to radiation has been neglected. 

(g) chemical production terms : 

Consider the chemical reactions 

S 3 

where A a represents the species s. Then the mass-production rate of 
species s is 

u>. = M, \k,AT) n^,r- - K tt (T) nix,]-*.' I (14) 

where the forward and backward production rates are expressed in the 
modified Arrhenius form as 

Kf k (T) = exp(Co,fc H — 7 jr- + C2,fclog e T) 

I<f k (T ) = exp(Do,fc + -pjr- + l°g e T ) 

where Co,*? . • • 5^2, it are given constants. 

Nondimensionalization 

All the physical quantities appearing in the preceding equations 
are nondimensionalized using their reference values . These reference 
values (dimensional quantities being denoted by *) are defined as : 

|r*U/ = V 

v ; e/ = c 

,* _ l r * |re 

l ref y * 

ref 

Pref Poo 


M‘„, = M" m 

n 


n» 

ref M* 


ref 


T . = Yik. 

K. t 

* T /* 2 

Pref Pref *ref 


* t /* 2 

e ref *ref 


Cl, r „ = Ktf 


Pref Poo 


K—f — Pref^ref 


Kef = 


W ref = 


Pref 

Pref 

P'r'fKf 


( 15 ) 


10 


and the other nondimensional parameter of interest is the Reynolds 
number which is defined as 


Re = 


pU f v~fV 

P re/ 


The equations Eq. 1, 2, and 3 along with the equations Eqs . 4-14 
provide a complete set of equations for the unlnown fields p 3 r p r \i r p r 
and T . These equations are discretized using the finite-volume ap- 
proach which is discussed next . 

Discretization 

(a) Preliminaries : 

The conservation equations, Eqs. 1, 2, 3 can be written as one sin- 
gle equation in nondimensional integral form as 


d_ 

dt 


f QdV + <b F ndA = / 
Jv(t) Jav(t) Mt) 


WdV 


(16) 


lv(t) JdV(t) 

where the vector of dependent conservative variables, Q f is 

T 

Q = {/M > P2. ■ ■ ■ > Pn. , P u 1 E] 


(17) 


F = F' - F 11 = 


Pi{ u-v) \ 


/ pWc\ \ 


( \ 

p 2 (u - v) 

1 

pVVc 2 

W = 

W 2 

Pn.( U-V) 

Re 

pVWc nt 



pu(u - v) + pi 
\ E( u — v) +pu / 


T 

\T • u — q/ 


0 

0 / 


(18) 


Let the n be the unit vector normal to the control surface. Fur 
ther, let F n = n-F be the normal component of the algebraic flux vector 
and let «„ = n-u, v n = n • v, < = u n - v n be respectively, the flow, the 


11 


control surface, and the relative velocity components normal to the 
control surface. The normal component of the inviscid flux vector is: 


T 

Fn(Q] n) = n • F* = {plU n , p 2 U n , • • • , Pn.u' n , P u n u + P n > ® U n + P u n} 

In many numerical algorithms, the above nonlinear flux vector is 
often linearized in time or space. This linearization leads to the 
following n e x n e inviscid Jacobian matrix 

u' n 6 rs — c r u n c r n 0 \ 


A* (17; n) = 


/ 

7(|u-u + ^)n-u„u un-7nu + u'„I jn 

V[ 7 (|u-u + 4> a )-H]u n Hn — 7« n u u' n +^u n / 
where S ra is the Kronecker delta with r, s = 1,2, . . . , n a and 

1 


(19) 


7 = 


, C Vf TM 

e “ 


( 20 ) 


MC V , 

U is the vector of primitive variables p a ,p,u,p, The eigenvalues of 

this inviscid Jacobian matrix are: 

A , = u_ i = 1,2,..., fi a -H Tiii 1 


An, + n<( — u n ®/ 

An, + n,i + l ”1” ®/ 

where the frozen speed of sound a/ is defined as 


( 21 ) 




( 22 ) 


Define the unit tangent vectors s, t such that n,s,t form a unit or 
thonormal spatial basis, i.e., 


n • s = n • t = s • t = 0 |n| — |s| — |t| — 1 
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also be written as A = RAL, where A is a 


The Jacobian matrix A can 
diagonal matrix of eigenvalues of A, 

A(*7;n) = diag {u' n ,u' n , . . • , u' n ,u' n - a f ,u' n + a f } T (23) 

and L and R are the left and right eigenvector matrices, respectively. 


R({7; n,s,t) = 


L(C7 ; n,s,t) = 


2a 


6 rs 

0 0 

C r 

C r N 

\ 

U 

aft a/s 

u — a/n 

u + a/n 

(24) 

|u • U - tps 

a/u • t a/u • s 

H — afu n 

H + afu n , 

/ 

/ 2a f &rs 

- 2 c r 7 (|u- u + r/> 3 ) 

2c r 7U 

- 2 c r 7 \ 



— 2a/u • t 

2a ft 

0 



— 2a/u • s 

2a fS 

0 

(25) 

7(2 u ' 

U + V’s) + a f u n 

— 711 — a/ 

n 7 


V 7(| u ' 

U + V’s) — a /“n 

-711 + a/ 

n 7 / 



Columns 2 and 3 of R and rows 2 and 3 of L can be dropped for the case of 
a one-dimensional flow. Column 3 of R and row 3 of L can be dropped for 
the case of a two-dimensional flow. 


(b) Spatial Discretization : 

In this research effort, the finite-volume method 7 is chosen to 
discretize the governing equations. In the finite-volume method the 
flow domain is divided into contiguous cells and the conservation prin 
ciples applied to each of these cells. This in effect replaces the 
surface integrals in Eq. 16 by the sum of integrals over the faces of 
the cells. For the case of three dimensional flow, the cells are hex- 
ahedra defined by specifying the vertices r j± i iJk± i 5 ,±i (whole indices 
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represent the cell centroid or cell center) and these vertices are 
used to compute the surface area vectors and cell volumes. In the in- 
terpretation of the volume integrals, the state of dependent vari- 
ables is assumed to be the value of Q at some average point in the cell, 
e.g., the cell center . For a hexahedral cell (see Fig. 1) centered at 
j,k,l, the semi-discrete version of Eq. 8 can be written as 


A nd 

+ E = v w»w 

a=l 


(26) 


m a represents either j , or k, or l according to whether a is 1, or 2, or 
3. In any case the index m a is always an integer. The circumflex above 
the flux and source vectors indicates that these are numerical quan- 
tities. These are consistent with the physical fluxes and source vec- 
tors. For example, if the numerical flux vector F a is defined as: 

F® ,1 = F a (Qm a +2i Qm a +li Qm a i Qm a — 1 > ^m a + 

771 a i 2 * 

then the consistency requirement is simply 

F a (Qm a 5 Qm a J Qm a i Qm a > S ma + l) — ) 

The required surface area vectors of the cell faces are computed 
using the formulas below 7 

S j±i,fc,i = \( r i±b k + l 2’ l -% “ ^ - r j±k,k-k>i-0 

±k t i = 2^ r i+b k± k’ l ~2 ~ r i-2>* ± 2> ,+ 5^ ® ( r i-b k± 2> l ~i r j+%,k±%,i+%) 

S j,fc,i±i = \( r j+b k +b l± h - T J-b k -b l± 0®( r i-b k +b l± ? ~ T j+b k -b l ±0 
and the cell volume is computed using the formula 7 

V,,M = i(r J+ i it+ 1,1+1 - 


(27) 


(28) 
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There are two ways of obtaining second and higher order representa 
tion of the inviscid numerical flux at a call face. In the first ap- 
proach, the cell centroid values at the neighboring cells are used 
to set up a Riemann problem at the face under consideration and this 
is called the MUSCL (Monotonic Upwind Schemes for Conservation Laws) 
approach. This approach is closer to the spirit of the finite-volume 
formulation. In the second approach, the physical fluxes at computed 
at the neighboring cell centers and then used in a weighted manner. 
This approach is called the non-MUSCL approach approach. Both these 
approaches are discussed below. 

( i ) Inviscid Fluxes - MUSCL representation 

In the MUSCL approach 15 the piecewise constant initial data of the 
Riemann problem in Godunov' s method 16 are replaced with piecewise lin 
ear inital data. The required right ( R ) and left (L) states of the 
Riemann problem at the face j + f a re defined as 



The parameter <f> determines the spatial order of accuracy. For exam 
pie, (f> = — 1 leads to a fully upwind scheme and 4> = \ leads to a third- 
order upwind— biased scheme. The left and right states as defined above 
lead to nonphysical oscillations at discontinuities. In order to elim- 
inate these oscillations we define the following intermediate vec- 
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tors 


Xl,m a + ± = (Qm a +2 — Qm a + 1) 

X2 ,m a + % = {Qm„ + \ ~ Qm a ) 

X3 ,m a + % ~ (Q™a ~ Qm a - 1) 

The elements of these intermediate vectors are then limited relative 
to each other. The limited vectors are given by 

Xi,m a +i = mimnod(xi ima+ i,6x 2 ,m 0 +i) 

X 2 ) m a + i = minmod(x2,m a + i> & X 3 ,m a + i) 
f 2 ,m 0 +i = minmod(x2, m< ,+i, 6 Xi,m a + i) 

%,m a + ± = minmod(x 3l m a + i J & X 2 ) m„ + i) 
where the limiting operator is defined by 


minmod(x, by) = sgn(x)max {0,min[|x|, 6ysgn(x)]} 

and 1 < b < (3 - <j>)/(l - <f>) and <f> ± 1 . In terms of these limited intermedi- 
ate vectors, the limited left and right states at a cell face are 


Qmc + i — ^ m * +1 


t+i - + 


1 — <t>_ 1 + 4> = 

^ Xl^a + i ' 4 X 2 ,m a + i 

1 — <f> = 1 + 4> — 

^ X3,m a + ^ ' 4 X2,m a + ^ 


Using these limited states in the first-order Roe scheme we obtain 
the higher-order inviscid numerical flux 

Ki+i = |{ s L + i - t Fi («“„ +i ) + F‘«?‘. +i )] - iA*-; +i iw2. +i ( 29 > 


where 


|A-: +i l = 5 L+t R» m „ +i |A” m „ +i |L:. +i 


(30) 


‘QTj 
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and the elements of the diagonal matrix |A| are 


l(A" +i ) ( |=tAP“„ + iW «= 1 , 2 ,...,"« 


The function t/>[z] is necessary to prevent entropy violating expansion 
shocks and defined as 11 


I (z 2 + e 2 )/2e |z| < e 


(31) 


and e is some small positive parameter. This function is required at 
points where the eigenvalue goes to zero. The superscript <* on L, R, 
and A indicates that these matrices are evaluated using the basis vec- 
tors n / s , t for the a-direction and the subscript m a + \ indicates that 
an averaged state at the cell face is used in the evaluation. This av 
eraged state is some symmetric average of the limited left and right 
states. Note that the spatial differencing stencil consists of five 
points . 

(ii ) Inviscid Fluxes - Non-MUSCL representation 

(1) Osher-Chakravarthy : 

The first-order accurate inviscid numerical flux is written 


as 


10 


cv 

TTO ot i 2 


= ^{s“ a+ i • [F < (Qm.+i)+F i (Q m J] -|A“'* + J(Qm Q+ i-Q m j[ (32) 


FO 


which is actually the average of two one-sided representations 


m a + 1 


F a ’ 1 


= S« Q + i • r(Q ma ) + + dQm a +l - Qm a ) 


FO 


m a + 


a i 2 


= S“ a + i • F'((5mc, + l) — A m ’^ii(Qm a +l Qm a ) 


FO 


a i 2 




(33) 
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where |A| is defined in Eq. 30 and the matrices A ± are 


A a,<± , = S a X .R“ ,,A a± , ,L“ , , 

m a + \ m a + \ m « + l m <* + 2 


The elements of the diagonal matrices A* are 




,n e 


where the function ip[z] is defined in Eq. 31 . 

Xn order to improve the spatial accuracy to second order the f ol 
lowing intermediate algebraic vectors (characteristic variables) 

are defined 3 _ _ 

X“m a +i = L l t> +iW n, «+ 2 “ Q"*a+l) 

X2,m a + i = L m a + ±(Q m ° + 1 ~ 

X3,m a + i = L m a + i(^« “ Qrn a -l) 

As before, to prevent nonphysical oscillations at shocks, the ele- 
ments of these intermediate vectors are limited relative to each other 
The limited vectors are 

X? >ma+ i =mininod(x“ ma+ i, 6 X 2 , ma+ i) 

X?, ma+ i = m inmod(x2, ma+ i^X3 )m « + i) 

+ i = minmod(x 2 , mtt + i > 6 x" ma + i) 

=minmod(x3 ima+ i,&X2,m a + i) 

In terms of these limited fluxes, the higher order inviscid flux can 
be written as 12 


rrOf,i 

' m Q + - 


m a + ^ 


l -4> 


ar 


FO 


R a ,xA , 

m a + i m a + 


iXi 
? 1 » 


T7loi + : 


A a,— V** i 

4 K m a + i A m£ , + i X 2,m Q + ^ 




1 + ^ 


A a,-f i 1 ^ r> & A a ’^ TT 01 

^m a + i^m 0 + 5 X3,m a + i + 4 m 0 + i m a + ^ 2 i m « + 5 


( 34 ) 
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The parameter <f> has the same meaning as in the MUSCL approach. Again, 
the matrices R, L, and A are evaluated at some averaged state at the 
cell face. This averaging could be done using Roe averaging proce- 
dure or even a simple arithmetic averaging procedure . Roe-averaging 
in the context of nonequilibrium is a little difficult since we can 
only satisfy either the jump condition or the jump condition but not 
both simultaneously. In the present work simple arithmetic averaging 
is used. Note that the spatial differencing stencil again consists of 
five points . 

(2) Harten: 

In this method of representing the numerical flux, the concept of 
a modified flux is used. In terms of this modified flux 11 , the second- 
order accurate inviscid numerical flux F ' , i is written as 

*71 ot "r j 


ma + 5 


HO 



[F i (g ma +i) + F i (Q ma )] 


+ S“ 


,,R fl 

1 2 


(Q 


m a + 1 



(35) 


The elements of the vector $ are 


Cti = + + 9ml) - + i ^ 


f 1 




where g l is the Ith element of the modified flux g and x l the Ith element 
of the characteristic variable x • As before the characteristic vari- 
able is defined as 


- I'“(Q nlo + i)(Qm 11 +l Qm a ) 


In the expression for the flux, the modified eigenvalue 7 1 is 

1 . f (0L a +l -?ma)/ a In a + i + 

= «L- 4 =0 
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and the modified flux is 


9 l m a =minmod(a' ma _i,a^ a+ i) 


Note that in this representation, we can achieve second order accu 
racy in both space and time. This scheme, however, is a little more 
dissipative than the other two schemes that have been outlined in the 
preceding sections . 

(iii) Viscous Fluxes 

The computation of the viscous part of the numerical flux involves 
the evaluation of the stress tensor and heat flux vector at the cell 
interfaces, i.e., we have to compute T moi± i , q m<I ±i • First consider 
the viscous stress tensor T ma+ i . This can be written as 




(Vu) m „ + i + (Vu)^ + j - j(V ■ u) m „ + il 


(36) 


In order to evaluate the terms indicated in the parentheses, we use 
the divergence theorem. 7 For an arbitrary field a (scalar or vector) 
we have 

/(V«a)^ = /a«ndA (37) 

where ® represents either a gradient, or a divergence, or a curl op- 
eration. Applying Eq. 36 to the gradient of the velocity field u on an 
auxiliary cell ABCD (seeFig. 2) surrounding the point indexed j + k, l 
(m a =j) yields 




u j+i,M S j+i,fc,J ” + u j+5.*+5> /S i+5,*+i,/ 


— Jk-i,/ + u j+b k ' l+ i^j+k’ k < l+ 2 
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The expression for the divergence of the velocity field can be ob- 
tained in a similar manner. Using these expressions, the final form 
of the viscous stress tensor at a cell face is 




a 

V 


j+b k > 1 


u i+ i,jfc,|S| +liJM - Uj t k,iS^ k t - S* j k l u jjkyl 


2 

— o ( u i+ 1 > fc > i ' S j+i,*,/ “ u i’ k < 1 ' 


u(\ 

?) f 

J V 



u i+ifc+i,i S i+i fc+i i 


- u i+b k -b‘ S Hb k - l 2’ 1 + U i+h k > l+ ± S< i+b k ’ l H U j+l k ’ l -?i+b k ’ l -h 

+ S T j+^k+bi Uj+ $’ k+ %' 1 ~ ST j+b k -b lUj+ %’ k ~%' 1 + S i+b k ’ l+ i u - ,+ 2- fc ’ /+ i 

2 

— S^ + i fc j_i u j+i,fc,i-i — 3 ( u >+|,*+^,i ‘ u j+b k ~b l ^j+b k ~b l 

+ u j +b k ’ l +i ' S i+i M+i ” u i+i. fc .'-| ' 


The first term in brackets represents the thin-layer contribution to 
the viscous flux and the second term in brackets represents the con- 
tributions due to the cross-derivatives. The geometrical quantities 
(volumes and area vectors) in the above expression are simple arith- 
metic averages. For the scalar and vector fields there are many dif- 
ferent ways of doing the averaging. Consider the velocity field u. 
The possible averaging procedures are 

( \ [uj+l,fc,/ + Uj+l,fc+l,f + Uj,fe+1,( + u j,k,l] 


U 


j+ i 


S 2 [ u i.M + u j+i,fc+i,i] 
k |[u i+ i,Jk ) / + u j,fc+l,f] 


The choice of averaging is dictated by diagonal dominance in an im- 
plicit scheme. Therefore, for an implicit scheme the last two choices 
improve diagonal dominance. The criteria for choosing the averaging 
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procedure are given in Ref. 12. The heat flux vector can be derived in 
a similar manner and the final form of this vector is given by 


l y 


j+l.M 




K 


y'i+\,k,n 


~^^J + ^ >M+ i ^ j+ \ ,M+ 5 ^3+ ^ i ^j+ \,k,l— \ 


The viscous numerical fluxes for the other faces and directions can 
be derived in a similar manner. 

(c) Temporal differencing : 

So far, the time-differencing part of the algorithm was not ad- 
dressed because we considered only the semi -discrete conservation 
law. If At represents the time steps, then the temporal term of Eq. 16 
is represented as 


4- [ QdV = [(1 + w)A n Q jlk ,i ~ uA n 1 Q i ,j k| /] 

dtJv(t) 

where A n Q jt k 1 i = Q]t)~Qlk,i and ^ controls the temporal accuracy of this 
representation. For example, if cj = 0 we have first-order temporal 
accuracy and if u> = | we have second-order temporal accuracy. 


Algorithm 


Using the concepts developed in the preceding sections, the numer- 
ical algorithm is written as a two-parameter family of implicit and 
explicit schemes. 

r n d 


A n Qj t k,i + 


At 


e 
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Cl? 


.— ip .. At (l- g ) V(P°.“ , ) - V, i, ,1V", , 


(38) 


(l+u>) 

where 0 = 0 defines an explicit scheme and 0^0 defines an implicit 
scheme. Further, second-order temporal accuracy is achieved only 
when 0 = uj + i . in the present study we concentrate only on implicit 
schemes. Before considering the actual solution procedure we assume 
that the grid is invariant in time, i.e.,v = 0. This also implies that 
the geometrical parameters such as the volume and surface area vec 


tors are invariant in time. 

The usual method of solving the nonlinear difference equation is 
through time linearization, i.e., the nonlinear equation is linearized 
about the known state Q n and the linear system of equations is solved 
using known techniques. This approach suffers from linearization er 
rors . In order to avoid linearization errors a Newton method is used 
As an analogy consider a simple nonlinear equation f(x ) = 0. The Newton 
method for this equation is 


f'(x k )(x k+1 - x k ) = -/(**) k = 1,2 , ... 


In the nonlinear difference equation represented by Eq. 38, Q n+1 is 
the unknown quantity. Let q p be an approximation to Q n+1 such that as 
p increases, q p approaches Q n+1 . Linearizing Eq. 38 about the known 
state q p we have 


At 


0 


n<i 


[ I+ viLattj ( S “s? **■'«« *J 


A p q = - Qlk,i)~ 


At 0 
Vj,k,l (1 + 


r rid 




O' — 1 


+ 


U) 


n — 1 


(1 +cj)' 


Qi,k,i 


2 3 


( 38 ) 


A t (1-0) 

Vj.fc,/ (1 +w) 


nd 

■a=l 


where A p q = q p+l - q p and 6 a () = () ma+ i - () m „ _i similarly the other terms . 
If q° = Q n and only one subiteration is carried out then we have the 


usual noniterative algorithm. The LHS of the algorithm can be sim 
plified by considering only a spatially first-order accurate scheme. 
Even so, when the subiterations converge, the RHS is satisfied to the 
desired accuracy. We also assume that eigenvalues and eigenvectors 
are independent of q p . Using the one-sided representation of the in- 
viscid numerical flux, Eq. 33, and making using of the fact that 2/ s/ = 


0 for a closed cell, we get 


i- Ate 


dW' 

dq. 


A p q + 


At e 
Vi,u( l+«) 


A*;+ (A%-A p qj -x)+ 

J 2 


A^; 4 (A p q j+1 - A %) + Al’+ h (A p q k - A p 9 *_i) + A*., (A p q k+1 - A p q k )+ 


A ( '.’\(A p qi — A p qi-i) + A^i(A p gi+i — A p qi) 

J o 2 


ip 


= RHS(Eq. 39) 


(40) 


The above equation is too large to solve on the computer and hence needs 
to be simplified. There are several options available to us . The first 


of these is to use a Gauss-Seidel relaxation in the predominant stream 
direction. 12 Let this direction be the j-direction. We then drop the 
j — 1 and j + 1 terms to obtain 


i — Ate 


dW 


+ 


9 ( t ,+ 

v J>M (i+u>r 



A p q j + 


At 


Al’\(A p q k - A p q k -i) + AJ’ i (A P 9 Jt + i - A p q k )+ 


Vj, k ,i C 1 + w ) 

A^l+ (A p <z, - A p qi-i) + A C '~(A p qi+i - A p q t ) 

J o J ' 2 


= RHS(Eq. 39) 


(41) 
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This equation is still too large to solve and is now impemented in an 
approximately factorized form. First define the operator A as 

dW At 9 


A* = 


t - (\t> + A*'") 

1 At ° dq + (1 + u;/ A M 


(42) 


Then the approximately factorized scheme is 

At 9 


(A .«)- 


A { + 


A< + (1 + u) 

At 9 


(A "’ + 6 b v + A *-6{) 


(A t'+S* + A <•-£*) A % = RHS(Eq. 39) (43) 


Vj,*,i(l+ W ) 

The operators S b amd 6 f represent the conventional backward and for- 
ward difference operators, respectively. Equation 43 is the final 
form of the algorithm which is implemented in the computer code. In 
its present form the algorithm is for three-dimensional flows. For a 
two-dimensional flow there is no approximate factorization and the 
algorithm reduces to the following line-relaxation algorithm 


A< + „ A ‘ - ~ ; (A"’ + ^ + W-%) 


A r gj = RHS(Eq. 39) 


(44) 


V>,», i(l+«) 

For a one-dimensional flow this further reduces to apoint-relaxatron 
scheme given by 

At' p A p qj = RHS(Eq. 39) (45) 


Note that for the two- and three-dimensional cases the operators rep- 
resent block-tridiagonal systems of equations where the order of the 
block matrices is n e . The one-dimensional case involves the inversion 
of a full matrix of order n t . 
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Boundary Conditions 


The discussion in this section is only applicable to external flows. 
For such a situation, the typical computational domain (s) (see Fig. 

3) are bounded by (i) an inf low boundary, (ii) an outflow boundary , 

(iii) a wall boundary, (iv) a freestream boundary, and (v) asymmetry 
boundary. There are additional special cases to consider such as (i) 
a geometrically singular boundary (axis) , and (ii) zonal or overlap 
boundaries . Each of these boundary conditions is discussed below : 
Inflow Boundary 

For the types of flows of interest, the inflow boundary is assumed 
to be in the freestream and Dirichlet boundary conditions are imple- 
mented implicitly. A special case of this boundary condition occurs 
in the three-dimensional flow around blunt bodies where the entire 
face collapses into a line . Such a boundary has to be treated with spe- 
cial care as it has been known to have caused problems in other three- 
dimensional calculations. This boundary is discussed later in this 
section. For a freestream inflow we have 

= Qoo VlT", k, l 


Outflow Boundary 

In the this research effort, the outflow is assumed to be super- 
sonic and hence a simple zeroth/first extrapolation of flow variables 
is done. This can be implemented implicitly. The extrapolation causes 


26 



some error in the subsonic layer close to the wall. However, past ex- 
perience has shown that this effect is negligible. This is expressed 

mathematically as 


Q 


n+l 

NJ,k,l 


= Q 


n+l 

NJ-l,k,l 


Vn, fc, l 


Wall Boundary 

The wall is assumed to noncatalytic, i.e., n • Vc a = 0 . Further the 

wall is assumed to be a isothermal with a fixed wall temperature T w 
which in principle can vary from point to point on the surface . No slip 
conditions u = 0 are also implemented at the wall. Assuming no-slip 
condition at the wall leaves us with only the pressure to be estimated 
at the wall . Assuming near boundary-layer behavior at the wall, this 
pressure is obtained by zeroth order extrapolation from the interior. 
Strictly speaking, the value of the wall pressure should be deter- 
mined by solving the normal momentum equation. 7 In the present case 
this is not done . 

The implementation of the wall boundary condition is quite involved. 
The layer of cells close to the wall have indices j',2,/ and cells in- 
dexed j,l,/ are fictitious cells (see Fig. 4) . The values of the flow 
variables in these cells are determined from the simple relations given 

below 


T”+; = 2 T. - T"+? 


In evaluating the viscous stresses, one sided differences are used at 
the wall . 
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Freestream Boundary 


Like the inflow boundary, at the freestream boundary Dirichlet con- 
ditions are used implicitly in the code. This far field boundary con- 
dition assumes that all flow discontinuities are contained well within 

this boundary . 

Qj#K,i = Q<» Vn ’ J ’ z 

Symmetry Boundary 

The three-dimensional geometry considered in this research are as- 
sumed to possess a plane of symmetry. Across this plane of symmetry, 
all flow quantities are reflected. These reflection boundary condi- 
tions are implemented implicitly and are expressed as 


T/n+l 

U j,k,NL 


— u n+1 

— u j,k,NL+l 


U — psi Pi Pi ^ 


w j,k,NL ~ W j,k,NL + 1 


U = p„p,p,T,E,u,v 

V = p„p,p,T,E 


Axis Boundary 13 

This boundary condition is perhaps the most difficult one to im- 
plement. It must be emphasized here that a geometric singularity does 
not imply that the flow variables are singular. There is some respite 
from this singularity in the finite-volume method because the cell 
centers done not coincide with the axis (see Fig. 5) . In the present 
analysis it is assumed that the grid is sufficiently smooth in the 
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neighborhood of the singularity. First we write the semidiscrete ver- 
sion of Eq. 1 for a cell indexed 1 ,k,l. However the state Q considered 
is assumed to be located at 5/6, k, l where 

5 1 

a similar conservation equation can be written for the call indexed 
Q 2kJ . Using this equation, Q 2 ,k,i can be eliminated from the conserva- 
tion equation for the cell indexed The variables required for 

the higher-order flux corrections are obtained from cells that are 
diametrically opposite to the cell being considered (see Fig. 5) . 

This is probably the single most important requirement. In this ef- 
fort we generate grids such that the singular axis is enclosed in a 
cylindrical tube extending two cell widths. This reduces the require- 
ment of interpolation to obtain the state vector at diametrically op- 
posite locations . 

Overlap Boundary 

Dirichlet boundary conditions are used at overlap boundaries . The 
values of Q from one grid are transferred to another through the over- 
lap region. If the overlapping regions do not exactly match then we 
have to use interpolation to transfer the values . The tranferred val 
ues are assumed to remain constant during the subiteration. 

Finally, it must be noted that low order boundary conditions are 
used in the implicit part of the algorithm. These boundary conditions 
are explicitly corrected to second-order accuracy at the end of every 
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subiteration, inasmuch as the left hand side goes to zero at conver- 
gence. We are also investigating the possibility of using character- 
istic boundary conditions. 

Code Development 

Before delving into the details of the new multidimensional code 
for reacting flows, the important aspect of grid generation has to be 
addressed. This is discussed next . 

(a) Grid Generation : 

The grid required by the new code is always three-dimensional ir 
respective of the dimensionality of the flow being computed. For a 
one-dimensional flow the grid is a tube of a unit square cross-section 
and for a two-dimensional flow the grid is a slab of unit depth. The 
reason for this is that a three-dimensional grid is necessary to com- 
pute the surface area vectors and cell volumes required by the finite- 
volume approach. In doing so one removes conditionality (which slows 
down the code) and the grid can be treated the same way for all flows . 
The price paid in this approach is increased storage for the one- and 
two-dimensional cases. In the present research, however, the very 
important aspect of grid generation is kept independent of the code, 
i.e.t the reacting flow code does not include any grid generation pack 
age. This makes the code more flexible and applicable to a variety of 

problems and faster at the same time. 

Now the code to integrate the equations is divided into three codes 
a pre-processor, an integrator, and a post-processor . Each of these 
codes is discussed below. 
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(b) Preprocessor : 


The basic unit in the integration of the equations is a sequence. 

A sequence can contain one or more finite-volume grids. Currently, 
these grids have to have a continuous overlap. For example, in the 
case of hypersonic flow past a wedge, one can have two sequences; the 
first sequence consisting of a single grid around the forebody and 
the second sequence containing two grids in the wake. Each sequence 
can have its own spatial and temporal accuracy. For the example of 
the wedge, the forebody sequence can have first-order temporal and 
second-order spatial accuracy, while the aftbody sequence can have 
second-order temporal and spatial accuracy. 

The function of the preprocessor is to set up the datasets and files 
necessary for the integration package. The preprocessor reads the 
flow conditions, the species set and the number of reactions. The pro 
gram then gets the species thermodynamic, transport and reaction data 
from the database. Currently, the data base contains the species 0 2 t 
N 2 , no, NO + , O, N, and e“ . The thermodynamic data are in the form 
cubic-spline coefficients . Using these data the code sets up the ref- 
erence conditions required for nondimensionalization . The number of 
sequences is then read in along with the pertinent accuarcy informa- 
tion. The coordinates of the grids that form the sequences are also 
read in. Based on these coordinates the surface area vectors and cell 
volumes are computed and written into files and the flow variables are 
initizalized at two time levels. The details of the reference con- 
ditions and file structure are written into the preprocessor output 
file which is read by the integration package. 
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(c) Integrator: 


This is the main code of the three and is a FORTRAN implementation 
of the algorithm shown in Eq. 43. It is in this code that the flow equa 
tions are actually integrated. The integration package reads in the 
reference conditions and file structure from the preprocessor file. 
Note that the integration package does not read in the grid coordi- 
nates but reads in the cell surface area vectors and volumes instead. 
This saves both computational time and storage since the area vectors 
and volumes are not computed over and over again and the grid coordi 
nates are not stored. This assumes that the grids used in the compu- 
tation are invariant in time. If at later date dynamic adpative grid- 
ding is desired then this code will have to be modified. 

It must be noted here that the boundary condition procedures are 
also not integral to the code. The code provides hooks for boundary 
conditions and it is up to the user to write the boundary condition 
procedures and routines. This was done in order to make the code flex- 
ible. Datasets in the code's native format and residual history are 
the only outputs from the integration package. These datasets and 
histories are analyzed by the postprocessor which is discussed next . 

(d) P ostprocessor 

The main function of the postprocessor is to analyze the f lowfield 
datasets computed by the integration package. The postprocessor com- 
putes the various fields like pressure, temperature, Mach number, 
etc. for purpose of graphical representation. The most important func 
tion of the postprocessor is the computation of the aerothermal loads. 
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The postprocessor computes the heat -transfer and aerodynamic coeffi- 
cients . The postprocessor can also be used for creating files for the 
purpose of computing radiation intensities and transport . 

Results 


(a) Model Problem 

In order to test the feasiblity of the numerical method outlined in 
the preceding sections, a model problem was solved. This model prob- 
lem is shown below in the integral form 

if udV+I ( U u+<l>i)dA=J ^dA+J u(u-\)(u-l)dV ( 46 ) 

dtj v(t) Lv(t) 2 Javtt)dx Jv(t) 2 

where w, ip, /i, « are constants . The second term on the left hand side of 
the equation represents the advection term. The first and second terms 
on the right hand side represent the viscous and source terms, re- 
spectively . The model problem is actually Burgers equation and in the 

differential form is 

Equation 4 6 represents many problems and of these a few representa- 
tive problems have been selected. These are 

Problem 1 : The first model problem solved was that of steady, in- 
viscid, non-reacting flow, i.e. , oj = 0, fi = k = 0, and ip = 1 . Mathemat 
ically stated, 

— + i-(—) = 0 

dt dx y 2 ' 

subject to the following initial and boundary conditions 

u(x,0) = l-2x 0 < x < 1, u(0,<) = M(M) = “I 0 <t < oo 
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The exact solution to the above problem is 

f 1 0 < x < 0.5 

u(x,t) - [_ x 0.5 < X < 1 

which means that there is a stationary shock located at x — 0.5. This 
case was computed using a CFL number of 1 x 10 6 . The calculations con- 
verged within four steps and the results are shown in Fig. 6. It is ev- 
ident that discontinuity that develops can be captured without os- 
cillations within two to three cells. Such crisp capture is due to the 
upwind nature of the scheme. 

Problem 2 : The second model problem solved was that of unsteady, 
inviscid, non-reacting flow. This was to check the time-accuracy (es- 
pecially important in the calculation of unsteady flows) of the algo- 
rithm. The problem statement is identical to the the one above except 
that the initial and boundary conditions are changed. Specifically, 
the initial and boundary conditions used are 


u 



0 < x < 0.3 
0.3 < x < 1 


tx(0,f) = l t >0 


The exact solution to the above problem is of a discontinuity travel 
ling to the right with a speed of stopped at t = 0.6, the discontinuity 
would be located at x = 0.6. The results of this calculation are shown 
in Fig. 7 . The figure clearly shows that the shock is tracked nearly 

perfectly . 

Problem 3 : 14 The third model problem solved was that of unsteady, 
inviscid, reacting flow. This was to check the influence of the source 
term on the calculation. The source term is somewhat contrived in the 
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sense that it does not represent the actual source term of a chemi- 
cally reacting system. The initial and boundary conditions used are 
the same as in the previous problem. The problem has been solved for 
various values of the source term strength k . The results of this cal- 
culation are shown in Figs. 8a-8f . Several important conclusions can 
be made from these computations. For small values of k , there is lit- 
tle or no influence on the solution. The source term tends to sharpen 
the discontinuity. When the source term reaches a critical value the 
discontinuity does not move very far from its initial position. It is 
unlikely that this situation will be encountered in the actual chem- 
ically reacting problem. A high source term strength translates to 
near equilibrium conditions. For such a situation it is not a good 
idea to use the nonequilibrium code to compute the equilibrium flow 
- one has to use a code dedicated to computing the equilibrium flow. 

Problem 4 : The final model problem solved was that of steady, vis 
cous, reacting and non-reacting flow. For the non-reacting case, the 
exact solution of the problem is known. The computed solution and the 
corresponding convergence history are shown in Figs. 9a-9b. Note that 
the reaction term again sharpens the gradient. 

These four model problems clearly demonstrate the feasilbilty and 
accuracy of the algorithm. The FORTRAN code developed for the actual 
chemically reacting system is currently being tested for a variety 
of problems. The results of these calculations will be presented at a 

later date . 
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Concluding Remarks 


The equations governing the reacting flow of a multicomponent gas 
have been derived. These equations are applicable to one-, two- and 
three-dimensional flows. The various modelling assumptions and ex- 
pressions have been detailed. A numerical scheme based on the finite- 
volume method has been developed to solve the governing equations . 

The various methods of representing the inviscid and viscous numer- 
ical fluxes have been detailed. The algorithm is spatially second 
order accurate (third order accuracy can be achieved in principle) 
and temporally second-order accurate. This temporal accuracy can 
be reduced to first order for cases where a steady state solution is 
required. The feasibility and accuracy of the algorithm hace been 
demonstrated for some simple one-dimensional model problems. For 
the actual problem of multi-dimensional chemically reacting flows, 
the algorithm has been implemented as a set of three computer codes 
a preprocessor, an integrator, and a postprocessor. The code is cu- 
urently being debugged and tested for simple geometrical shapes in 
two dimensions. Once the code has been validated for the simpler cases, 
it will be used to compute the complete unsteady reacting flow field 
around the Aeroassist Flight Experiment (AFE) body shape. 
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Figure 1. Typical finite volume cell showing cell corners 
(open circles) and cell center (solid circle) 



Figure 2 . Secondary grid for transport calculations 
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Figure 3. Typical multiple grid system for external flows 
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Figure 8a. Solution with no source term. 







Figure 8d. Solution with source strength=10 . 0 










Figure 9a. Solution of the viscous, non-reacting problem 
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